Disclaimer

This work is very preliminary as I get back into the coding swing of things. Data wrangling and figure generation will be done via R, but the rest of the project will be done using good ol’ microsoft products. This is just an entry point into data crunching and should by no means be considered a final product.

3/30/2023 update, corrected a minor error in methodology.

Steamboat COOP

Coop Data

NOAA/NWS Cooperative Observer Network

Two stations near steamboat, according to this site. 057936 02 SBTC2 UNITED STATES CO ROUTT +7 STEAMBOAT SPRINGS 40 30 17 -106 51 58 6636 057942 02 SSPC2 UNITED STATES CO ROUTT +7 STEAMBOAT SPRINGS 1 W 40 29 00 -106 51 00 6700

Downloading from the Iowa Environmental Mesonet. NOAA package not working, maybe due to user error.

Data Cleaning

#coop_steamboat_clean <- read.csv("C:/Users/bsteen/OneDrive - #Colostate/Documents/MS_Research/research-COOP/data_clean/coop_steamboat_clean.csv", header = TRUE)

Figure check

# from home
coop_steamboat_clean <- read.csv("C:/Users/13074/Documents/ESS580/thesis_project/research-COOP/data_clean/coop_steamboat_clean.csv", header = TRUE)

coop_steamboat_clean$Date <- ymd(coop_steamboat_clean$Date)

str(coop_steamboat_clean)
## 'data.frame':    44740 obs. of  13 variables:
##  $ ï..station  : chr  "CO7936" "CO7936" "CO7936" "CO7936" ...
##  $ station_name: chr  "STEAMBOAT SPRINGS" "STEAMBOAT SPRINGS" "STEAMBOAT SPRINGS" "STEAMBOAT SPRINGS" ...
##  $ lat         : num  40.5 40.5 40.5 40.5 40.5 ...
##  $ lon         : num  -107 -107 -107 -107 -107 ...
##  $ day         : chr  "1/1/1900" "1/2/1900" "1/3/1900" "1/4/1900" ...
##  $ doy         : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ max_T_c     : num  -5.6 -0.6 5 7.8 3.3 4.4 -1.7 0.6 0.6 0.6 ...
##  $ min_T_c     : num  -21.1 -17.8 -6.1 -5.6 -10.6 -9.4 -18.9 -15.6 -15 -18.9 ...
##  $ avg_T_c     : num  -13.35 -9.2 -0.55 1.1 -3.65 ...
##  $ Date        : Date, format: "1900-01-01" "1900-01-02" ...
##  $ waterYear   : int  1900 1900 1900 1900 1900 1900 1900 1900 1900 1900 ...
##  $ daymonth    : chr  "01-01" "02-01" "03-01" "04-01" ...
##  $ waterDay    : int  93 94 95 96 97 98 99 100 101 102 ...
ggplot(coop_steamboat_clean, aes(x = Date, y = avg_T_c)) +
  geom_line() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

#Check for outliers....

#dygraph

COOP_temp_xts <- xts(coop_steamboat_clean$avg_T_c, order.by = coop_steamboat_clean$Date)

dygraph(COOP_temp_xts) %>%
  dyAxis("y", label = "Daily temperature (°C)") 

Detrending Data

#SF figured out the yearly average by water year

#average water year temperature

COOP_yearly_wy_aver <- coop_steamboat_clean %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

COOP_daily_wy_aver <- COOP_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

COOP_daily_wy_aver <- COOP_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(COOP_daily_wy_aver$aver_day_temp))

#str(daily_wy_aver)
# try to show all years as means. 
COOP_daily_wy_aver2 <- COOP_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
COOP_daily_wy_aver2$date_temp <- signif(COOP_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(COOP_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard Deviation

To figure out the standard deviation for each year, I want the “residual” for each daily value.

The standard deviation will be the daily residual minus the mean of the residuals by water year, summed and squared, then divided by the number of observations minus one. The square root of the resulting value of which is thus the standard deviation for the water year.

Determining residuals

COOP_standard_dev <- COOP_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(COOP_standard_dev$residual)
## [1] -9.652271e-17

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

COOP_standard_dev_all <- COOP_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

COOP_standard_dev_all <- COOP_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

COOP_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1900 3.387043
1901 4.068619
1902 3.568775
1903 3.527189
1904 3.696915
1905 3.402598
1906 3.167953
1907 3.996705
1908 3.453319
1909 4.792216
1910 4.499426
1911 4.143812
1912 3.935198
1913 4.109351
1914 4.003080
1915 3.687804
1916 3.984756
1917 4.499360
1918 3.987734
1919 3.751334
1920 4.167794
1921 4.301257
1922 4.401713
1923 3.695604
1924 3.607483
1925 3.947664
1926 3.400954
1927 4.236689
1928 4.342036
1929 4.295103
1930 4.208991
1931 3.797143
1932 4.291793
1933 4.295393
1934 3.508106
1935 3.647796
1936 3.725678
1937 3.945438
1938 3.655086
1939 4.076417
1940 3.720026
1941 3.833928
1942 4.269962
1943 4.039773
1944 3.637570
1945 3.828605
1946 3.936963
1947 3.833509
1948 4.127381
1949 3.747441
1950 4.092183
1951 4.305097
1952 4.137060
1953 4.104637
1954 3.493113
1955 4.223100
1956 4.473028
1957 4.227564
1958 3.727149
1959 3.861405
1960 3.846644
1961 3.447334
1962 4.549574
1963 4.254752
1964 4.116082
1965 4.272005
1966 3.505496
1967 3.952601
1968 3.786388
1969 4.124864
1970 4.031999
1971 4.382322
1972 4.300176
1973 4.211447
1974 3.657256
1975 4.026061
1976 3.425550
1977 3.371104
1978 3.745753
1979 4.228923
1980 3.970125
1981 3.800191
1982 4.062310
1983 3.823264
1984 3.937347
1985 4.034121
1986 3.904996
1987 3.411860
1988 3.964910
1989 4.381722
1990 3.725159
1991 4.245862
1992 4.030424
1993 4.086199
1994 3.757401
1995 4.037364
1996 3.935450
1997 4.092659
1998 3.733417
1999 4.194017
2000 3.619707
2001 4.252626
2002 4.292146
2003 3.756071
2004 3.940469
2005 3.678190
2006 4.100780
2007 4.276960
2008 3.721782
2009 3.681365
2010 3.902448
2011 4.177220
2012 3.455744
2013 4.447458
2014 3.878538
2015 4.109502
2016 3.583526
2017 4.422096
2018 3.662640
2019 4.031328
2020 4.427930
2021 3.851773
2022 4.567149
ggplot(COOP_standard_dev_all, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Update 5/4/2023

write.csv(COOP_standard_dev_all,"C:/Users/13074/Documents/ESS580/thesis_project/research-COOP/data_clean/COOP_standard_dev_all.csv", row.names = FALSE)
# THIS IS A NIGHTMARE (but it worked).

auto_cor_SD <-COOP_standard_dev_all %>% 
  mutate(Decade = paste0(waterYear  %/% 10 * 10))

acf_1900 <- auto_cor_SD %>%
  filter(Decade == "1900")
acf_1910 <- auto_cor_SD %>%
  filter(Decade == "1910")
acf_1920 <- auto_cor_SD %>%
  filter(Decade == "1920")
acf_1930 <- auto_cor_SD %>%
  filter(Decade == "1930")
acf_1940 <- auto_cor_SD %>%
  filter(Decade == "1940")
acf_1950 <- auto_cor_SD %>%
  filter(Decade == "1950")
acf_1960 <- auto_cor_SD %>%
  filter(Decade == "1960")
acf_1970 <- auto_cor_SD %>%
  filter(Decade == "1970")
acf_1980 <- auto_cor_SD %>%
  filter(Decade == "1980")
acf_1990 <- auto_cor_SD %>%
  filter(Decade == "1990")
acf_2000 <- auto_cor_SD %>%
  filter(Decade == "2000")
acf_2010 <- auto_cor_SD %>%
  filter(Decade == "2010" | Decade == "2020")

acf(acf_1900$sd_2, lag=5, pl=FALSE)
## 
## Autocorrelations of series 'acf_1900$sd_2', by lag
## 
##      0      1      2      3      4      5 
##  1.000 -0.239  0.202 -0.160 -0.207 -0.052
acf(acf_1910$sd_2, lag=5, pl=FALSE)
## 
## Autocorrelations of series 'acf_1910$sd_2', by lag
## 
##      0      1      2      3      4      5 
##  1.000  0.037 -0.533  0.129  0.143 -0.320
acf(acf_1920$sd_2, lag=5, pl=FALSE)
## 
## Autocorrelations of series 'acf_1920$sd_2', by lag
## 
##      0      1      2      3      4      5 
##  1.000  0.231 -0.047 -0.218 -0.470 -0.284
acf(acf_1930$sd_2, lag=5, pl=FALSE)
## 
## Autocorrelations of series 'acf_1930$sd_2', by lag
## 
##      0      1      2      3      4      5 
##  1.000  0.018 -0.087  0.018 -0.113 -0.277
acf(acf_1940$sd_2, lag=5, pl=FALSE)
## 
## Autocorrelations of series 'acf_1940$sd_2', by lag
## 
##      0      1      2      3      4      5 
##  1.000 -0.095 -0.486 -0.097  0.019  0.166
acf(acf_1950$sd_2, lag=5, pl=FALSE)
## 
## Autocorrelations of series 'acf_1950$sd_2', by lag
## 
##      0      1      2      3      4      5 
##  1.000  0.076 -0.525 -0.442  0.286  0.283
acf(acf_1960$sd_2, lag=5, pl=FALSE)
## 
## Autocorrelations of series 'acf_1960$sd_2', by lag
## 
##      0      1      2      3      4      5 
##  1.000 -0.135 -0.049 -0.187 -0.402  0.151
acf(acf_1970$sd_2, lag=5, pl=FALSE)
## 
## Autocorrelations of series 'acf_1970$sd_2', by lag
## 
##      0      1      2      3      4      5 
##  1.000  0.414  0.088 -0.178 -0.206 -0.464
acf(acf_1980$sd_2, lag=5, pl=FALSE)
## 
## Autocorrelations of series 'acf_1980$sd_2', by lag
## 
##      0      1      2      3      4      5 
##  1.000 -0.054 -0.526 -0.001  0.162 -0.116
acf(acf_1990$sd_2, lag=5, pl=FALSE)
## 
## Autocorrelations of series 'acf_1990$sd_2', by lag
## 
##      0      1      2      3      4      5 
##  1.000 -0.537  0.190 -0.422  0.462 -0.293
acf(acf_2000$sd_2, lag=5, pl=FALSE)
## 
## Autocorrelations of series 'acf_2000$sd_2', by lag
## 
##      0      1      2      3      4      5 
##  1.000 -0.083 -0.484 -0.073 -0.032  0.430
acf(acf_2010$sd_2, lag=5, pl=FALSE)
## 
## Autocorrelations of series 'acf_2010$sd_2', by lag
## 
##      0      1      2      3      4      5 
##  1.000 -0.610  0.383 -0.113  0.043 -0.087
#put the output into an excel file, reading it back into R. #failingbutwinning


steamboat_autocorrelation <- read_excel("data_clean/steamboat_autocorrelation.xlsx")
#View(steamboat_autocorrelation)

steamboat_autocorrelation <- steamboat_autocorrelation %>% 
  filter(Decade != "2010_2022")

ggplot(steamboat_autocorrelation, aes(x = Decade))+
  geom_point(aes(y = lag_1, color = "lag_1"), size =7)+
  geom_point(aes(y = lag_2, color = "lag_2"), size =5)+
  geom_point(aes(y = lag_3, color = "lag_3"), size =4)+
  theme_few()+
  scale_colour_manual(name = "lags", values=c("lag_1"="red", "lag_2" = "purple", "lag_3" ="blue"))+#, "noaa_morrisey" = "red", "oyler" = "purple"))+#, "oyler_noaa_conus" = "grey"))+
  xlab("Decade")+
  ylab("Autocorrelation")

coop_steamboat_clean_annual_t <- coop_steamboat_clean %>% 
  group_by(waterYear) %>% 
  mutate(mean_temp_annual = mean(avg_T_c)) %>%  
  distinct(waterYear, .keep_all = TRUE)

ggplot(coop_steamboat_clean_annual_t, aes(x = waterYear, y = mean_temp_annual)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Mean annual temperature (°C)') + 
  xlab('Water Year')

Checking a random year from the timeseries.

COOP_standard_dev_87 <- COOP_standard_dev %>% 
  filter(waterYear == 1987) %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
           mutate(sd_1 = residual-resid_mean)

COOP_standard_dev_87 <- COOP_standard_dev_87 %>%
  group_by(waterYear) %>%
  mutate(sd_2 = (((sum((sd_1)^2))/((sum(tabulate(COOP_standard_dev_87$waterDay)))-1)))^(0.5))

head(COOP_standard_dev_87$sd_2, 1)
## [1] 3.41186

Looks good.

Summer temperature standard deviation

COOP_standard_dev_all_summer <- COOP_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

COOP_standard_dev_all_summer <- COOP_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

COOP_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1900 2.343534
1901 2.052822
1902 2.370343
1903 2.038227
1904 1.629967
1905 1.889724
1906 1.768985
1907 1.296196
1908 1.748968
1909 2.170115
1910 2.000957
1911 2.208233
1912 2.288062
1913 2.526328
1914 1.713976
1915 2.080550
1916 2.122182
1917 1.921234
1918 2.714718
1919 2.404358
1920 2.379640
1921 2.272216
1922 1.883791
1923 1.584776
1924 2.169737
1925 2.168214
1926 2.394504
1927 2.673759
1928 2.450568
1929 2.395691
1930 1.470897
1931 2.339435
1932 2.157538
1933 1.947806
1934 1.711129
1935 1.957617
1936 2.086720
1937 2.119000
1938 1.939949
1939 2.255753
1940 2.115228
1941 1.983932
1942 1.681608
1943 2.104935
1944 1.844256
1945 2.463376
1946 1.699111
1947 2.293489
1948 2.131599
1949 1.329732
1950 1.680547
1951 2.123582
1952 1.787384
1953 2.011197
1954 2.419949
1955 2.080150
1956 1.824454
1957 2.189850
1958 1.994329
1959 2.140958
1960 1.620575
1961 1.606953
1962 2.063112
1963 1.526357
1964 2.460590
1965 1.479728
1966 1.769717
1967 1.429726
1968 2.292879
1969 2.257943
1970 1.902575
1971 1.787828
1972 1.928531
1973 2.586347
1974 2.239603
1975 1.953987
1976 2.229307
1977 1.959354
1978 1.810978
1979 1.875752
1980 1.908943
1981 2.143438
1982 2.335733
1983 2.053450
1984 1.967233
1985 2.174145
1986 1.863319
1987 1.904726
1988 1.986655
1989 2.285310
1990 2.617198
1991 1.793621
1992 2.324413
1993 2.552794
1994 2.116211
1995 2.749449
1996 1.750358
1997 2.340247
1998 2.430688
1999 1.866542
2000 2.215133
2001 2.618700
2002 2.528510
2003 2.415446
2004 2.541391
2005 2.222956
2006 2.165075
2007 2.405073
2008 2.596269
2009 1.977769
2010 2.115651
2011 1.943490
2012 2.072429
2013 1.943989
2014 2.253113
2015 2.337471
2016 1.947347
2017 1.914768
2018 2.316994
2019 2.483465
2020 2.575885
2021 2.482696
2022 3.147948
ggplot(COOP_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Jun-Aug standard deviation for water years 1900-2021

Mann-Kendall & Sen’s Slope

Summer standard deviations.

sd_mk_summer <- mk.test(COOP_standard_dev_all_summer$sd_2)
print(sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_all_summer$sd_2
## z = 2.2211, n = 123, p-value = 0.02635
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 1.017000e+03 2.092503e+05 1.355458e-01
sd_sens_summer <- sens.slope(COOP_standard_dev_all_summer$sd_2)
print(sd_sens_summer)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_all_summer$sd_2
## z = 2.2211, n = 123, p-value = 0.02635
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.0002669607 0.0035814520
## sample estimates:
## Sen's slope 
## 0.001931784

Winter temperature standard deviation

COOP_standard_dev_all_winter <- COOP_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# COOP_standard_dev_all_winter <- COOP_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


COOP_standard_dev_all_winter <- COOP_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

COOP_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1900 4.032437
1901 5.152863
1902 4.580459
1903 4.296317
1904 4.277880
1905 4.516199
1906 4.042376
1907 3.615871
1908 4.196627
1909 6.027292
1910 6.176796
1911 5.179486
1912 5.285784
1913 5.418424
1914 5.124939
1915 4.492648
1916 5.320355
1917 5.960599
1918 5.006565
1919 4.710630
1920 5.211839
1921 5.465862
1922 6.082555
1923 4.710294
1924 4.409977
1925 4.966388
1926 3.926105
1927 5.161938
1928 5.548213
1929 5.336150
1930 5.268585
1931 4.785408
1932 5.872585
1933 5.506158
1934 4.296671
1935 4.823304
1936 4.478443
1937 5.261662
1938 4.633655
1939 5.180047
1940 4.921106
1941 5.032277
1942 5.812571
1943 5.213421
1944 4.709201
1945 4.595392
1946 5.257146
1947 4.751165
1948 5.401913
1949 5.149603
1950 5.149961
1951 5.986132
1952 5.734292
1953 5.455903
1954 4.601240
1955 5.585458
1956 6.282696
1957 5.781879
1958 5.015481
1959 4.849098
1960 5.105714
1961 4.087986
1962 6.156861
1963 6.117983
1964 5.350773
1965 5.941975
1966 4.716248
1967 5.124600
1968 4.841351
1969 5.489272
1970 5.297339
1971 5.753360
1972 6.016482
1973 5.445182
1974 4.815183
1975 5.371131
1976 4.575222
1977 4.206960
1978 4.424031
1979 5.922070
1980 5.425155
1981 4.268802
1982 5.399768
1983 4.560962
1984 5.047083
1985 5.158547
1986 5.208834
1987 4.194948
1988 5.013068
1989 5.757657
1990 4.458001
1991 5.520284
1992 4.848221
1993 5.208729
1994 4.739457
1995 5.039713
1996 4.769257
1997 4.751369
1998 4.682027
1999 5.320874
2000 4.075872
2001 4.969522
2002 5.145239
2003 4.611630
2004 4.854507
2005 4.674537
2006 5.038415
2007 5.375765
2008 4.262245
2009 4.469792
2010 4.613645
2011 5.463104
2012 4.168804
2013 5.434713
2014 4.645119
2015 5.265263
2016 4.279881
2017 5.543005
2018 4.391600
2019 4.616723
2020 4.637958
2021 4.022003
2022 5.033060
ggplot(COOP_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Nov-Mar standard deviation for water years 1900-2021

Mann-Kendall & Sen’s Slope

Winter standard deviations.

sd_mk_winter <- mk.test(COOP_standard_dev_all_winter$sd_2)
print(sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_all_winter$sd_2
## z = -0.99685, n = 123, p-value = 0.3188
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
## -4.570000e+02  2.092503e+05 -6.090897e-02
sd_sens_winter <- sens.slope(COOP_standard_dev_all_winter$sd_2)
print(sd_sens_winter)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_all_winter$sd_2
## z = -0.99685, n = 123, p-value = 0.3188
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.004753212  0.001474537
## sample estimates:
##  Sen's slope 
## -0.001523523

Smaller time increments might be useful….

1987 - 2021

In comparison to the Tower SNOTEL site:

# I did this wrong- need to compare the 1987- 2021 time series without the mean of the 1900 - 2021 time series influencing the residuals.  

#average water year temperature

COOP_yearly_wy_aver_87_21 <- coop_steamboat_clean %>%
  filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))

#Average temperature by day for all water years:

COOP_daily_wy_aver_87_21 <- COOP_yearly_wy_aver_87_21 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

COOP_daily_wy_aver_87_21 <- COOP_daily_wy_aver_87_21 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(COOP_daily_wy_aver_87_21$aver_day_temp)) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
  

COOP_standard_dev_87_21 <- COOP_daily_wy_aver_87_21 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

COOP_standard_dev_87_21 <- COOP_standard_dev_87_21 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)


COOP_standard_dev_87_21 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 3.424576
1988 3.931157
1989 4.301842
1990 3.583093
1991 4.176727
1992 3.937969
1993 3.973356
1994 3.691415
1995 3.964424
1996 3.939454
1997 4.033360
1998 3.667026
1999 4.157142
2000 3.663586
2001 4.116828
2002 4.183094
2003 3.693382
2004 3.871310
2005 3.644129
2006 4.118883
2007 4.051204
2008 3.733130
2009 3.738246
2010 3.926954
2011 4.218237
2012 3.297969
2013 4.352540
2014 3.875402
2015 4.077057
2016 3.566496
2017 4.354637
2018 3.658564
2019 3.989471
2020 4.365474
2021 3.782126
ggplot(COOP_standard_dev_87_21, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

SD for water years 1987 - 2021

1987 - 2021 Mann-Kendall & Sen’s Slope

sd_mk_87_21 <- mk.test(COOP_standard_dev_87_21$sd_2)
print(sd_mk_87_21)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_87_21$sd_2
## z = 0.68167, n = 35, p-value = 0.4954
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 4.900000e+01 4.958333e+03 8.235294e-02
sd_sens_87_21 <- sens.slope(COOP_standard_dev_87_21$sd_2)
print(sd_sens_87_21)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_87_21$sd_2
## z = 0.68167, n = 35, p-value = 0.4954
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.007011434  0.013956979
## sample estimates:
## Sen's slope 
## 0.003231668

Summer and Winter 87-21

Summer temperature standard deviation

# using the 1987- 2021 data frame

COOP_standard_dev_all_87_21_summer <- COOP_daily_wy_aver_87_21 %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual))) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

COOP_standard_dev_all_87_21_summer <- COOP_standard_dev_all_87_21_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

COOP_standard_dev_all_87_21_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 1.935143
1988 1.993794
1989 2.230940
1990 2.544550
1991 1.775141
1992 2.332864
1993 2.568171
1994 2.178839
1995 2.786218
1996 1.789394
1997 2.339336
1998 2.382172
1999 1.828553
2000 2.201437
2001 2.627660
2002 2.511709
2003 2.405932
2004 2.511997
2005 2.219459
2006 2.134931
2007 2.358026
2008 2.534519
2009 1.916377
2010 2.109093
2011 1.851202
2012 1.991655
2013 1.970984
2014 2.214680
2015 2.302278
2016 1.998966
2017 1.877157
2018 2.294063
2019 2.407073
2020 2.540686
2021 2.469644
ggplot(COOP_standard_dev_all_87_21_summer, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Summer season SD for water years 1987 - 2021

Mann-Kendall & Sen’s Slope

Summer 87-21 standard deviations.

sd_mk_summer <- mk.test(COOP_standard_dev_all_87_21_summer$sd_2)
print(sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_all_87_21_summer$sd_2
## z = 0.17042, n = 35, p-value = 0.8647
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 1.300000e+01 4.958333e+03 2.184874e-02
sd_sens_summer <- sens.slope(COOP_standard_dev_all_summer$sd_2)
print(sd_sens_summer)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_all_summer$sd_2
## z = 2.2211, n = 123, p-value = 0.02635
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.0002669607 0.0035814520
## sample estimates:
## Sen's slope 
## 0.001931784

Winter temperature standard deviation

# using the 1987- 2021 data frame

COOP_standard_dev_all_87_21_winter <- COOP_daily_wy_aver_87_21 %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual))) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

COOP_standard_dev_all_87_21_winter <- COOP_standard_dev_all_87_21_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

COOP_standard_dev_all_87_21_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 4.169279
1988 4.979472
1989 5.615506
1990 4.256383
1991 5.449026
1992 4.557310
1993 4.986382
1994 4.662775
1995 4.923980
1996 4.826996
1997 4.671113
1998 4.633265
1999 5.189225
2000 4.159790
2001 4.845489
2002 5.107031
2003 4.524114
2004 4.696476
2005 4.530915
2006 5.149943
2007 5.022096
2008 4.291589
2009 4.508637
2010 4.696043
2011 5.614273
2012 3.965244
2013 5.404382
2014 4.665509
2015 5.102774
2016 4.276957
2017 5.371856
2018 4.382774
2019 4.707970
2020 4.582198
2021 4.031273
ggplot(COOP_standard_dev_all_87_21_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Winter season SD for water years 1987 - 2021

Mann-Kendall & Sen’s Slope

1987 - 2021 winter Mann-Kendall & Sen’s Slope

sd_mk_87_21_winter <- mk.test(COOP_standard_dev_all_87_21_winter$sd_2)
print(sd_mk_87_21_winter)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_all_87_21_winter$sd_2
## z = -0.82368, n = 35, p-value = 0.4101
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -59.00000000 4958.33333333   -0.09915966
sd_sens_87_21_winter <- sens.slope(COOP_standard_dev_all_87_21_winter$sd_2)
print(sd_sens_87_21_winter)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_all_87_21_winter$sd_2
## z = -0.82368, n = 35, p-value = 0.4101
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.024144341  0.009118474
## sample estimates:
##  Sen's slope 
## -0.007927524

Neither the 1987 - 2021 as a whole, nor the seasonal time series demonstrate trends that are statistically significant.

1900-1915 minimum and maximum temperatures

Jan 7, 1913 had the lowest recorded temperature. This just looks at the max and min temperatures for that early part of the record.

early_coop <- coop_steamboat_clean %>% 
  filter(waterYear >= 1900 & waterYear <= 1915)

COOP_temp_min_xts <- xts(early_coop$min_T_c, order.by = early_coop$Date)

dygraph(COOP_temp_min_xts) %>%
  dyAxis("y", label = "Daily min temperature (°C)") 

Daily minimum temperatures for water years 1900-1915

COOP_temp_max_xts <- xts(early_coop$max_T_c, order.by = early_coop$Date)

dygraph(COOP_temp_max_xts) %>%
  dyAxis("y", label = "Daily max temperature (°C)") 

Daily minimum temperatures for water years 1900-1915

By groups of decades

1900-1930

#average water year temperature

COOP_yearly_wy_aver_00_30 <- coop_steamboat_clean %>%
  filter(waterYear >= 1900 & waterYear <= 1930) %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))

#Average temperature by day for all water years:

COOP_daily_wy_aver_00_30 <- COOP_yearly_wy_aver_00_30 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

COOP_daily_wy_aver_00_30 <- COOP_daily_wy_aver_00_30 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(COOP_daily_wy_aver_00_30$aver_day_temp)) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
  

COOP_standard_dev_00_30 <- COOP_daily_wy_aver_00_30 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

COOP_standard_dev_00_30 <- COOP_standard_dev_00_30 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)


COOP_standard_dev_00_30 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1900 3.201760
1901 3.990374
1902 3.495081
1903 3.498700
1904 3.566867
1905 3.350260
1906 3.103802
1907 3.966733
1908 3.464956
1909 4.752127
1910 4.430479
1911 4.027024
1912 4.005247
1913 4.093351
1914 3.911415
1915 3.677228
1916 3.900792
1917 4.528786
1918 3.859437
1919 3.921495
1920 4.094066
1921 4.162064
1922 4.342919
1923 3.735643
1924 3.683247
1925 3.807423
1926 3.383065
1927 3.993153
1928 4.145520
1929 4.180380
1930 4.206931
ggplot(COOP_standard_dev_00_30, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

SD for water years 1900 - 1930

1900 - 1930 Mann-Kendall & Sen’s Slope

sd_mk_00_30 <- mk.test(COOP_standard_dev_00_30$sd_2)
print(sd_mk_00_30)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_00_30$sd_2
## z = 2.1415, n = 31, p-value = 0.03223
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  127.0000000 3461.6666667    0.2731183
sd_sens_00_30 <- sens.slope(COOP_standard_dev_00_30$sd_2)
print(sd_sens_00_30)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_00_30$sd_2
## z = 2.1415, n = 31, p-value = 0.03223
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.001562166 0.032696841
## sample estimates:
## Sen's slope 
##  0.01672178

Summer and Winter 00-30

Summer temperature standard deviation

# using the 1900- 1930 data frame

COOP_standard_dev_all_00_30_summer <- COOP_daily_wy_aver_00_30 %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual))) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

COOP_standard_dev_all_00_30_summer <- COOP_standard_dev_all_00_30_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

COOP_standard_dev_all_00_30_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1900 2.263739
1901 2.090944
1902 2.264615
1903 2.122756
1904 1.608957
1905 1.891343
1906 1.751173
1907 1.447252
1908 1.762553
1909 2.266312
1910 2.006374
1911 2.164735
1912 2.296291
1913 2.560634
1914 1.738455
1915 2.069813
1916 2.063861
1917 1.898723
1918 2.568409
1919 2.358741
1920 2.317614
1921 2.145945
1922 1.821744
1923 1.557476
1924 2.135836
1925 2.059975
1926 2.348773
1927 2.590941
1928 2.453216
1929 2.387445
1930 1.490211
ggplot(COOP_standard_dev_all_00_30_summer, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Summer season SD for water years 1900 - 1930

Mann-Kendall & Sen’s Slope

Summer 00-30 standard deviations.

sd_mk_00_30_summer <- mk.test(COOP_standard_dev_all_00_30_summer$sd_2)
print(sd_mk_00_30_summer)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_all_00_30_summer$sd_2
## z = 1.3257, n = 31, p-value = 0.1849
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   79.0000000 3461.6666667    0.1698925
sd_sens_00_30_summer <- sens.slope(COOP_standard_dev_all_00_30_summer$sd_2)
print(sd_sens_00_30_summer)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_all_00_30_summer$sd_2
## z = 1.3257, n = 31, p-value = 0.1849
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.005609245  0.021399935
## sample estimates:
## Sen's slope 
## 0.007379067

Winter temperature standard deviation

# using the 00-30 data frame

COOP_standard_dev_all_00_30_winter <- COOP_daily_wy_aver_00_30 %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual))) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

COOP_standard_dev_all_00_30_winter <- COOP_standard_dev_all_00_30_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

COOP_standard_dev_all_00_30_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1900 3.723162
1901 5.065432
1902 4.492115
1903 4.195997
1904 4.199892
1905 4.449922
1906 3.941312
1907 3.778798
1908 4.266280
1909 5.919534
1910 6.051690
1911 5.107299
1912 5.454531
1913 5.346476
1914 4.990005
1915 4.512612
1916 5.220994
1917 5.903666
1918 4.886774
1919 4.984953
1920 5.066982
1921 5.216569
1922 5.999659
1923 4.691835
1924 4.563404
1925 4.706124
1926 3.993937
1927 4.787616
1928 5.195486
1929 5.033618
1930 5.213335
ggplot(COOP_standard_dev_all_00_30_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Winter season SD for water years 1900 - 1930

Mann-Kendall & Sen’s Slope

1900 - 1930 winter Mann-Kendall & Sen’s Slope

sd_mk_00_30_winter <- mk.test(COOP_standard_dev_all_00_30_winter$sd_2)
print(sd_mk_00_30_winter)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_all_00_30_winter$sd_2
## z = 1.5297, n = 31, p-value = 0.1261
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   91.0000000 3461.6666667    0.1956989
sd_sens_00_30_winter <- sens.slope(COOP_standard_dev_all_00_30_winter$sd_2)
print(sd_sens_00_30_winter)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_all_00_30_winter$sd_2
## z = 1.5297, n = 31, p-value = 0.1261
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.00447963  0.04646028
## sample estimates:
## Sen's slope 
##   0.0224029

1931-1960

#average water year temperature

COOP_yearly_wy_aver_31_60 <- coop_steamboat_clean %>%
  filter(waterYear >= 1931 & waterYear <= 1960) %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))

#Average temperature by day for all water years:

COOP_daily_wy_aver_31_60 <- COOP_yearly_wy_aver_31_60 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

COOP_daily_wy_aver_31_60 <- COOP_daily_wy_aver_31_60 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(COOP_daily_wy_aver_31_60$aver_day_temp)) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
  

COOP_standard_dev_31_60 <- COOP_daily_wy_aver_31_60 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

COOP_standard_dev_31_60 <- COOP_standard_dev_31_60 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)


COOP_standard_dev_31_60 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1931 3.760890
1932 4.193367
1933 4.319782
1934 3.529914
1935 3.567877
1936 3.702520
1937 3.878129
1938 3.610735
1939 3.971156
1940 3.672531
1941 3.800795
1942 4.099793
1943 3.912735
1944 3.570191
1945 3.735226
1946 3.906807
1947 3.798591
1948 4.062537
1949 3.668341
1950 4.116322
1951 4.225246
1952 4.056279
1953 4.107455
1954 3.577426
1955 4.177570
1956 4.385685
1957 4.170128
1958 3.652907
1959 3.800659
1960 3.774516
ggplot(COOP_standard_dev_31_60, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

SD for water years 1931 - 1960

1931 - 1960 Mann-Kendall & Sen’s Slope

sd_mk_31_60 <- mk.test(COOP_standard_dev_31_60$sd_2)
print(sd_mk_31_60)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_31_60$sd_2
## z = 1.3202, n = 30, p-value = 0.1868
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   75.0000000 3141.6666667    0.1724138
sd_sens_31_60 <- sens.slope(COOP_standard_dev_31_60$sd_2)
print(sd_sens_31_60)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_31_60$sd_2
## z = 1.3202, n = 30, p-value = 0.1868
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.003720908  0.021202461
## sample estimates:
## Sen's slope 
## 0.008983505

Summer and Winter 31-60

Summer temperature standard deviation

# using the 1931- 1960 data frame

COOP_standard_dev_all_31_60_summer <- COOP_daily_wy_aver_31_60 %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual))) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

COOP_standard_dev_all_31_60_summer <- COOP_standard_dev_all_31_60_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

COOP_standard_dev_all_31_60_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1931 2.233385
1932 2.110648
1933 1.930662
1934 1.641400
1935 1.888024
1936 2.023637
1937 2.102115
1938 1.945271
1939 2.301620
1940 2.008671
1941 1.896990
1942 1.696607
1943 2.023789
1944 1.894943
1945 2.383864
1946 1.708431
1947 2.286007
1948 2.194445
1949 1.312206
1950 1.738952
1951 1.999415
1952 1.768064
1953 1.950363
1954 2.430058
1955 1.958513
1956 1.838985
1957 2.258270
1958 2.032796
1959 2.024376
1960 1.568366
ggplot(COOP_standard_dev_all_31_60_summer, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Summer season SD for water years 1931 - 1960

Mann-Kendall & Sen’s Slope

Summer 31-60 standard deviations.

sd_mk_31_60_summer <- mk.test(COOP_standard_dev_all_31_60_summer$sd_2)
print(sd_mk_31_60_summer)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_all_31_60_summer$sd_2
## z = -0.24977, n = 30, p-value = 0.8028
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -15.00000000 3141.66666667   -0.03448276
sd_sens_31_60_summer <- sens.slope(COOP_standard_dev_all_31_60_summer$sd_2)
print(sd_sens_31_60_summer)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_all_31_60_summer$sd_2
## z = -0.24977, n = 30, p-value = 0.8028
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.014148582  0.007630278
## sample estimates:
##  Sen's slope 
## -0.002408317

Winter temperature standard deviation

# using the 31-60 data frame

COOP_standard_dev_all_31_60_winter <- COOP_daily_wy_aver_31_60 %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual))) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

COOP_standard_dev_all_31_60_winter <- COOP_standard_dev_all_31_60_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

COOP_standard_dev_all_31_60_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1931 4.751342
1932 5.790312
1933 5.562620
1934 4.284776
1935 4.748514
1936 4.474569
1937 5.162471
1938 4.573951
1939 5.084584
1940 4.793605
1941 4.987704
1942 5.525195
1943 5.081059
1944 4.556901
1945 4.504257
1946 5.245315
1947 4.679191
1948 5.324723
1949 5.087298
1950 5.122150
1951 5.929128
1952 5.695226
1953 5.489698
1954 4.703118
1955 5.587329
1956 6.162403
1957 5.712359
1958 4.909039
1959 4.761019
1960 4.948515
ggplot(COOP_standard_dev_all_31_60_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Winter season SD for water years 1931 - 1960

Mann-Kendall & Sen’s Slope

1931 - 1960 winter Mann-Kendall & Sen’s Slope

sd_mk_31_60_winter <- mk.test(COOP_standard_dev_all_31_60_winter$sd_2)
print(sd_mk_31_60_winter)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_all_31_60_winter$sd_2
## z = 1.6414, n = 30, p-value = 0.1007
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   93.0000000 3141.6666667    0.2137931
sd_sens_31_60_winter <- sens.slope(COOP_standard_dev_all_31_60_winter$sd_2)
print(sd_sens_31_60_winter)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_all_31_60_winter$sd_2
## z = 1.6414, n = 30, p-value = 0.1007
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.003118133  0.042830343
## sample estimates:
## Sen's slope 
##   0.0197381

1961-1990

#average water year temperature

COOP_yearly_wy_aver_61_90 <- coop_steamboat_clean %>%
  filter(waterYear >= 1961 & waterYear <= 1990) %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))

#Average temperature by day for all water years:

COOP_daily_wy_aver_61_90 <- COOP_yearly_wy_aver_61_90 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

COOP_daily_wy_aver_61_90 <- COOP_daily_wy_aver_61_90 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(COOP_daily_wy_aver_61_90$aver_day_temp)) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
  

COOP_standard_dev_61_90 <- COOP_daily_wy_aver_61_90 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

COOP_standard_dev_61_90 <- COOP_standard_dev_61_90 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)


COOP_standard_dev_61_90 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1961 3.356405
1962 4.532893
1963 4.290470
1964 4.043002
1965 4.287876
1966 3.519792
1967 3.881385
1968 3.715950
1969 4.106720
1970 3.940231
1971 4.283287
1972 4.232475
1973 4.157176
1974 3.552421
1975 3.917275
1976 3.416435
1977 3.309884
1978 3.721035
1979 4.065413
1980 3.968734
1981 3.755229
1982 3.899905
1983 3.798763
1984 3.970423
1985 3.972664
1986 3.782944
1987 3.502178
1988 3.887587
1989 4.262324
1990 3.756356
ggplot(COOP_standard_dev_61_90, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

SD for water years 1961 - 1990

1961 - 1990 Mann-Kendall & Sen’s Slope

sd_mk_61_90 <- mk.test(COOP_standard_dev_61_90$sd_2)
print(sd_mk_61_90)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_61_90$sd_2
## z = -1.1775, n = 30, p-value = 0.239
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##           S        varS         tau 
##  -67.000000 3141.666667   -0.154023
sd_sens_61_90 <- sens.slope(COOP_standard_dev_61_90$sd_2)
print(sd_sens_61_90)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_61_90$sd_2
## z = -1.1775, n = 30, p-value = 0.239
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.022821843  0.005520893
## sample estimates:
##  Sen's slope 
## -0.009830486

Summer and Winter 61-90

Summer temperature standard deviation

# using the 1961- 1990 data frame

COOP_standard_dev_all_61_90_summer <- COOP_daily_wy_aver_61_90 %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual))) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

COOP_standard_dev_all_61_90_summer <- COOP_standard_dev_all_61_90_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

COOP_standard_dev_all_61_90_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1961 1.603809
1962 2.004051
1963 1.463219
1964 2.423854
1965 1.530696
1966 1.759707
1967 1.376443
1968 2.256148
1969 2.240113
1970 1.813610
1971 1.826700
1972 1.927688
1973 2.561721
1974 2.295467
1975 1.892478
1976 2.192216
1977 1.940471
1978 1.799323
1979 1.858812
1980 1.875592
1981 2.094430
1982 2.355018
1983 1.999217
1984 1.991654
1985 2.149504
1986 1.868181
1987 1.954196
1988 1.986891
1989 2.275018
1990 2.559757
ggplot(COOP_standard_dev_all_61_90_summer, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Summer season SD for water years 1961 - 1990

Mann-Kendall & Sen’s Slope

Summer 61-90 standard deviations.

sd_mk_61_90_summer <- mk.test(COOP_standard_dev_all_61_90_summer$sd_2)
print(sd_mk_61_90_summer)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_all_61_90_summer$sd_2
## z = 2.0339, n = 30, p-value = 0.04196
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  115.0000000 3141.6666667    0.2643678
sd_sens_61_90_summer <- sens.slope(COOP_standard_dev_all_61_90_summer$sd_2)
print(sd_sens_61_90_summer)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_all_61_90_summer$sd_2
## z = 2.0339, n = 30, p-value = 0.04196
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.001372517 0.024724573
## sample estimates:
## Sen's slope 
##  0.01391241

Winter temperature standard deviation

# using the 61-90 data frame

COOP_standard_dev_all_61_90_winter <- COOP_daily_wy_aver_61_90 %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual))) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

COOP_standard_dev_all_61_90_winter <- COOP_standard_dev_all_61_90_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

COOP_standard_dev_all_61_90_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1961 3.915725
1962 6.144011
1963 6.187684
1964 5.258583
1965 5.976798
1966 4.779479
1967 5.036343
1968 4.757923
1969 5.468078
1970 5.116522
1971 5.610921
1972 5.910732
1973 5.367749
1974 4.623904
1975 5.213719
1976 4.574714
1977 4.053940
1978 4.422553
1979 5.671879
1980 5.426113
1981 4.231795
1982 5.115570
1983 4.505216
1984 5.153583
1985 5.080134
1986 5.027172
1987 4.379475
1988 4.862919
1989 5.556552
1990 4.526207
ggplot(COOP_standard_dev_all_61_90_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Winter season SD for water years 1961 - 1990

Mann-Kendall & Sen’s Slope

1961 - 1990 winter Mann-Kendall & Sen’s Slope

sd_mk_61_90_winter <- mk.test(COOP_standard_dev_all_61_90_winter$sd_2)
print(sd_mk_61_90_winter)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_all_61_90_winter$sd_2
## z = -1.8911, n = 30, p-value = 0.0586
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##           S        varS         tau 
## -107.000000 3141.666667   -0.245977
sd_sens_61_90_winter <- sens.slope(COOP_standard_dev_all_61_90_winter$sd_2)
print(sd_sens_61_90_winter)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_all_61_90_winter$sd_2
## z = -1.8911, n = 30, p-value = 0.0586
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.052990614  0.002647224
## sample estimates:
## Sen's slope 
## -0.02459531

1991-2020

#average water year temperature

COOP_yearly_wy_aver_91_20 <- coop_steamboat_clean %>%
  filter(waterYear >= 1991 & waterYear <= 2020) %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))

#Average temperature by day for all water years:

COOP_daily_wy_aver_91_20 <- COOP_yearly_wy_aver_91_20 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

COOP_daily_wy_aver_91_20 <- COOP_daily_wy_aver_91_20 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(COOP_daily_wy_aver_91_20$aver_day_temp)) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
  

COOP_standard_dev_91_20 <- COOP_daily_wy_aver_91_20 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

COOP_standard_dev_91_20 <- COOP_standard_dev_91_20 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)


COOP_standard_dev_91_20 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1991 4.145178
1992 3.960420
1993 3.952914
1994 3.687460
1995 3.931778
1996 3.931379
1997 4.023023
1998 3.624635
1999 4.094488
2000 3.636966
2001 4.116507
2002 4.208402
2003 3.707167
2004 3.870138
2005 3.641583
2006 4.115911
2007 4.058783
2008 3.757877
2009 3.718922
2010 3.931119
2011 4.226286
2012 3.267981
2013 4.357008
2014 3.840112
2015 4.062662
2016 3.616594
2017 4.330801
2018 3.629519
2019 3.966571
2020 4.346426
ggplot(COOP_standard_dev_91_20, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

SD for water years 1991 - 2020

1991 - 2020 Mann-Kendall & Sen’s Slope

sd_mk_91_20 <- mk.test(COOP_standard_dev_91_20$sd_2)
print(sd_mk_91_20)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_91_20$sd_2
## z = 0.28546, n = 30, p-value = 0.7753
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 1.700000e+01 3.141667e+03 3.908046e-02
sd_sens_91_20 <- sens.slope(COOP_standard_dev_91_20$sd_2)
print(sd_sens_91_20)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_91_20$sd_2
## z = 0.28546, n = 30, p-value = 0.7753
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.009547233  0.014518748
## sample estimates:
## Sen's slope 
## 0.001987054

Summer and Winter 91-20

Summer temperature standard deviation

# using the 1991- 2020 data frame

COOP_standard_dev_all_91_20_summer <- COOP_daily_wy_aver_91_20 %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual))) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

COOP_standard_dev_all_91_20_summer <- COOP_standard_dev_all_91_20_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

COOP_standard_dev_all_91_20_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1991 1.789832
1992 2.359212
1993 2.532099
1994 2.202599
1995 2.726217
1996 1.829452
1997 2.321800
1998 2.308894
1999 1.830620
2000 2.179309
2001 2.606471
2002 2.519656
2003 2.373394
2004 2.524740
2005 2.217524
2006 2.152131
2007 2.361085
2008 2.517364
2009 1.889837
2010 2.127676
2011 1.800363
2012 2.030146
2013 1.961039
2014 2.194757
2015 2.357020
2016 2.075895
2017 1.871864
2018 2.322724
2019 2.363316
2020 2.457244
ggplot(COOP_standard_dev_all_91_20_summer, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Summer season SD for water years 1991 - 2020

Mann-Kendall & Sen’s Slope

Summer 91-20 standard deviations.

sd_mk_91_20_summer <- mk.test(COOP_standard_dev_all_91_20_summer$sd_2)
print(sd_mk_91_20_summer)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_all_91_20_summer$sd_2
## z = -0.6066, n = 30, p-value = 0.5441
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -35.00000000 3141.66666667   -0.08045977
sd_sens_91_20_summer <- sens.slope(COOP_standard_dev_all_91_20_summer$sd_2)
print(sd_sens_91_20_summer)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_all_91_20_summer$sd_2
## z = -0.6066, n = 30, p-value = 0.5441
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.016065617  0.007924092
## sample estimates:
##  Sen's slope 
## -0.003377995

Winter temperature standard deviation

# using the 91-20 data frame

COOP_standard_dev_all_91_20_winter <- COOP_daily_wy_aver_91_20 %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual))) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

COOP_standard_dev_all_91_20_winter <- COOP_standard_dev_all_91_20_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

COOP_standard_dev_all_91_20_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1991 5.396681
1992 4.580670
1993 4.981946
1994 4.620652
1995 4.894775
1996 4.826410
1997 4.651361
1998 4.557633
1999 5.099055
2000 4.117501
2001 4.813144
2002 5.120805
2003 4.571242
2004 4.673569
2005 4.517890
2006 5.137815
2007 5.033545
2008 4.361708
2009 4.495860
2010 4.700523
2011 5.632693
2012 3.900892
2013 5.441360
2014 4.616701
2015 5.074498
2016 4.353704
2017 5.348487
2018 4.344638
2019 4.663561
2020 4.607316
#This was trying to help Nick Chohan out with a figure issue. Fail!

coef(lm(waterYear ~ sd_2, data= COOP_standard_dev_all_91_20_winter))
## (Intercept)        sd_2 
## 2014.582247   -1.903548
#str(COOP_standard_dev_all_91_20_winter)
 
ggplot(COOP_standard_dev_all_91_20_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_abline(xintercept = 1990, slope = 0.01333333, size =1) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') +
  xlab('Water year')

Winter season SD for water years 1991 - 2020

Mann-Kendall & Sen’s Slope

1991 - 2020 winter Mann-Kendall & Sen’s Slope

sd_mk_91_20_winter <- mk.test(COOP_standard_dev_all_91_20_winter$sd_2)
print(sd_mk_91_20_winter)
## 
##  Mann-Kendall trend test
## 
## data:  COOP_standard_dev_all_91_20_winter$sd_2
## z = -0.78501, n = 30, p-value = 0.4325
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -45.0000000 3141.6666667   -0.1034483
sd_sens_91_20_winter <- sens.slope(COOP_standard_dev_all_91_20_winter$sd_2)
print(sd_sens_91_20_winter)
## 
##  Sen's slope
## 
## data:  COOP_standard_dev_all_91_20_winter$sd_2
## z = -0.78501, n = 30, p-value = 0.4325
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02269601  0.01097902
## sample estimates:
##  Sen's slope 
## -0.005677511

Comparing Climate Normal means

(from 356:)

1991-2022

COOP_yearly_wy_aver_1991_2022_CN <- coop_steamboat_clean %>%
  filter(waterYear >= 1991 & waterYear <= 2022) %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))

#Average temperature by day for all water years:

COOP_daily_wy_aver_1991_2022_CN <- COOP_yearly_wy_aver_1991_2022_CN %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

COOP_daily_wy_aver_1991_2022_CN <- COOP_daily_wy_aver_1991_2022_CN %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(COOP_daily_wy_aver_1991_2022_CN$aver_day_temp)) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
COOP_daily_wy_aver_1991_2022_CN2 <- COOP_daily_wy_aver_1991_2022_CN %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))%>% 
  select(waterDay, date_temp) %>% 
  distinct(waterDay, .keep_all = TRUE)
  
COOP_daily_wy_aver_1991_2022_CN2$date_temp <- signif(COOP_daily_wy_aver_1991_2022_CN2$date_temp,3) #reduce the sig figs

ggplot(COOP_daily_wy_aver_1991_2022_CN2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  theme_few()

1961-1990

COOP_yearly_wy_aver_1961_1990_CN <- coop_steamboat_clean %>%
  filter(waterYear >= 1961 & waterYear <= 1990) %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))

#Average temperature by day for all water years:

COOP_daily_wy_aver_1961_1990_CN <- COOP_yearly_wy_aver_1961_1990_CN %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

COOP_daily_wy_aver_1961_1990_CN <- COOP_daily_wy_aver_1961_1990_CN %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(COOP_daily_wy_aver_1961_1990_CN$aver_day_temp)) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
COOP_daily_wy_aver_1961_1990_CN2 <- COOP_daily_wy_aver_1961_1990_CN %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))%>% 
  select(waterDay, date_temp) %>% 
  distinct(waterDay, .keep_all = TRUE)
  
COOP_daily_wy_aver_1961_1990_CN2$date_temp <- signif(COOP_daily_wy_aver_1961_1990_CN2$date_temp,3) #reduce the sig figs

ggplot(COOP_daily_wy_aver_1961_1990_CN2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  theme_few()

1931-1960

COOP_yearly_wy_aver_1931_1960_CN <- coop_steamboat_clean %>%
  filter(waterYear >= 1931 & waterYear <= 1960) %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))

#Average temperature by day for all water years:

COOP_daily_wy_aver_1931_1960_CN <- COOP_yearly_wy_aver_1931_1960_CN %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

COOP_daily_wy_aver_1931_1960_CN <- COOP_daily_wy_aver_1931_1960_CN %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(COOP_daily_wy_aver_1931_1960_CN$aver_day_temp)) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
COOP_daily_wy_aver_1931_1960_CN2 <- COOP_daily_wy_aver_1931_1960_CN %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))%>% 
  select(waterDay, date_temp)  %>% 
  distinct(waterDay, .keep_all = TRUE)
  
COOP_daily_wy_aver_1931_1960_CN2$date_temp <- signif(COOP_daily_wy_aver_1931_1960_CN2$date_temp,3) #reduce the sig figs

ggplot(COOP_daily_wy_aver_1931_1960_CN2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  theme_few()

1900-1930

COOP_yearly_wy_aver_1901_1930_CN <- coop_steamboat_clean %>%
  filter(waterYear >= 1901 & waterYear <= 1930) %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))

#Average temperature by day for all water years:

COOP_daily_wy_aver_1901_1930_CN <- COOP_yearly_wy_aver_1901_1930_CN %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

COOP_daily_wy_aver_1901_1930_CN <- COOP_daily_wy_aver_1901_1930_CN %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(COOP_daily_wy_aver_1901_1930_CN$aver_day_temp)) %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
COOP_daily_wy_aver_1901_1930_CN2 <- COOP_daily_wy_aver_1901_1930_CN %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))%>% 
  select(waterDay, date_temp) %>% 
  distinct(waterDay, .keep_all = TRUE)


COOP_daily_wy_aver_1901_1930_CN2$date_temp <- signif(COOP_daily_wy_aver_1901_1930_CN2$date_temp,3) #reduce the sig figs

ggplot(COOP_daily_wy_aver_1901_1930_CN2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  theme_few()

climate_norm <- left_join(COOP_daily_wy_aver_1991_2022_CN2, COOP_daily_wy_aver_1961_1990_CN2, by= 'waterDay') 

climate_norm <- left_join(climate_norm, COOP_daily_wy_aver_1931_1960_CN2, by= 'waterDay')

climate_norm <- left_join(climate_norm, COOP_daily_wy_aver_1901_1930_CN2, by= 'waterDay')

ggplot(climate_norm, aes(x=waterDay)) +  
         #scale_size_manual(values=c("FALSE"=5,"TRUE"=8)) +
  #scale_color_manual(values=c("FALSE"='#CC0000',"TRUE"='black')) +
geom_line(aes(y = date_temp.x, color = "1991-2022"), size =1)+
  geom_line(aes(y = date_temp.y, color = "1961-1990"), size =1)+
  geom_line(aes(y = date_temp.x.x, color = "1931-1960"), size =1)+
  geom_line(aes(y = date_temp.y.y, color= "1901-1930"), size =1)+
  theme_few()+
  scale_colour_manual(name = "climate normals", values=c("1991-2022" = "orange", "1961-1990" = "darkgreen", "1931-1960" = "purple", "1901-1930" = "red"))+
  xlab("Day of Year")+
  ylab("Average temperature")